Tidal Heating of Exomoons in Resonance and Implications for Detection

Mean motion resonance of celestial bodies in orbit around a primary has been studied in depth, with numerous applications ranging from solar system objects (Peale 1976) to specific exoplanetary systems (Correia et al. 2018). Here, we mostly follow the mathematical formulism set forth by Petrovich et al. (2013).

Our system consists of a planet as the primary, a test-particle moon in orbit around the planet, and an outer perturbing moon in a 2:1 resonance with the inner moon. Given this resonance, we have a relation between the period of the moons,

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{2}}{{P}_{1}}=\displaystyle \frac{2}{1},\end{eqnarray} \tag{ 1 }$$

where P₁ and P₂ are the orbital periods of the inner and outer moon, respectively. The Hamiltonian of the system is derived for this restricted three-body problem and can be written as

$$\begin{eqnarray}&&K=-3{\rm{\Delta }}R+{R}^{2}-2\sqrt{2R}\ \cos (r)\end{eqnarray} \tag{ 2 }$$

(see Petrovich et al. 2013 for full details of the calculation). Here, Δ is the resonance distance,

$$\begin{eqnarray}&&{\rm{\Delta }}=0.4{\left(\displaystyle \frac{{m}_{p}}{{m}_{2}}\right)}^{2/3}\left[2{\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)}^{5/6}-\displaystyle \frac{6+5{e}_{1}^{2}}{6+9{e}_{1}^{2}}{\left(\displaystyle \frac{{a}_{2}}{{a}_{1}}\right)}^{2/3}\right],\end{eqnarray} \tag{ 3 }$$

for planet mass m_p , inner moon with semimajor axis a₁ and small eccentricity e₁, and perturbing moon of mass m₂ and semimajor axis a₂. R is the canonical momentum,

$$\begin{eqnarray}&&R=0.9{\left(\displaystyle \frac{{m}_{p}}{{m}_{2}}\right)}^{2/3}{\left(\displaystyle \frac{{a}_{2}}{{a}_{1}}\right)}^{2/3}\left(\displaystyle \frac{6{e}_{1}^{2}}{6+7{e}_{1}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

and r is the canonical coordinate,

$$\begin{eqnarray}&&r=-2{\lambda }_{2}+{\lambda }_{1}+{\omega }_{1},\end{eqnarray} \tag{ 5 }$$

where λ₁ and λ₂ are the mean longitude of the inner and outer moon, respectively, and ω₁ is the longitude of periapsis of the inner moon. For constant Hamiltonian K, the semimajor axis and eccentricity of the inner moon will oscillate, or librate, about certain equilibrium values. These values are obtained by minimizing the Hamiltonian, which translates to solving the cubic equation

$$\begin{eqnarray}&&{x}^{3}-3{\rm{\Delta }}x-2=0,\end{eqnarray} \tag{ 6 }$$

where x is the canonical momentum in Cartesian coordinates,

$$\begin{eqnarray}&&x=\sqrt{2R}\ \cos (r).\end{eqnarray} \tag{ 7 }$$

Together with y, the canonical coordinate in Cartesian coordinates,

$$\begin{eqnarray}&&y=\sqrt{2R}\ \sin (r),\end{eqnarray} \tag{ 8 }$$

these coordinates track the orbital oscillations for constant K, as shown in Figure 1. In this case, arbitrary values of 0.5 and 1 are chosen for K and Δ, respectively, with Δ representing the shape of the moon trajectory in the constant Hamiltonian phase space. Because these values do not correspond to a minimum of K, there is libration of the orbital parameters, as can be seen by the bean-shaped trajectory. The libration period can be estimated to be (Murray & Dermott 1999)

$$\begin{eqnarray}&&{t}_{\mathrm{lib}}=0.4{P}_{1}{\left(\displaystyle \frac{{a}_{2}}{{a}_{1}}\right)}^{2/3}{\left(\displaystyle \frac{{m}_{p}}{{m}_{2}}\right)}^{2/3}\left(1+\displaystyle \frac{5}{6}{e}_{1}^{2}\right),\end{eqnarray} \tag{ 9 }$$

where P₁ is the orbital period of the inner moon. This shows that, as the planet mass dominates the perturbing moon or we move to larger MMR ratios, there is a slower change in the orbital values of the system.

Zoom In Zoom Out Reset image size

Thus, the closer to true resonance, the smaller the Hamiltonian of the system, and as it approaches a minimum the semimajor axis and eccentricity vary less and less (Murray & Dermott 1999). Conversely, the Hamiltonian can be increased to obtain a greater variation in orbit, while the moons are still considered to be in resonance. Nevertheless, for any degree of resonance we find a nonzero eccentricity is inevitable, which motivates our later calculations of tidal heating. In the example above, if the semimajor axis and eccentricity were chosen such that K was minimized, we would end up in one of the fixed points (−1, 0) or (2, 0) shown in red, and there would be no oscillation in these parameters (no libration). In other words, the system would be in perfect resonance.

A well-known example of a system in close to perfect resonance is Jupiter orbited by Io and Europa with orbital periods of 1.77 days and 3.55 days, respectively. Although the mass of Io is not negligible compared to the mass of Europa (a ratio of about 1/2) so that the restricted three-body problem is not directly applicable, we show that the formulas described so far can still be used to approximate the trajectory (Figure 2). For this system, Δ = −2.245 and K = −0.276, and because we are close to the Hamiltonian minimum the phase curve is simply a small circle around the fixed point. Thus, although we are not strictly in the restricted three-body regime, the method can still be applied for a test particle that is small compared to the central body.

Zoom In Zoom Out Reset image size

The equations above apply to two moons in the specific 2:1 resonance configuration, but they can be easily extended to apply to other first-order resonances with minor changes (see Petrovich et al. 2013 for details). Thus, because 2:1 mean motion resonance is simplest and most common to exoplanet systems (Lissauer et al. 2011), we choose to focus on this configuration. We will briefly discuss the effect other orders of resonance have on tidal heating in the Section 3.1.

Due to the complicated nature of even a three-body system under gravitational influence, we utilize the code REBOUND (Rein & Liu 2012), which is an N-body integrator that solves the motions of particles under the force of gravity. We first consider the zero eccentricity case, and take a generic three-body system consisting of a planet and two moons. The mass ratio of planet to outer moon is 300, and the moons are placed in a 2:1 resonance. Integrating forward in time using REBOUND, we obtain the trajectory shown in the left panel of Figure 3. Here, Δ = 0 and K = 0 because of the initial zero eccentricity, the location of this initial configuration is the origin of the x–y plot. The right panel of the figure shows the theoretical calculation of this system (as described in Section 2.1), which is similar, but not identical, since the simulated phase curve is smaller. This is due to the fact that REBOUND takes into account the center of mass of the system that the bodies orbit. The analytic model, by contrast, does not include these effects but has the center of mass located exactly at the origin of the coordinate system, the center of the planet. To test this explanation, we ran additional simulations for this system but with increasing mass ratio m_p /m₂. The result was a trend toward better agreement between theory and simulation for the higher mass ratio, confirming that the discrepancy is due the location of the center of mass. Nonetheless, the theory and numerics agree well even for a low mass ratio, demonstrating the robustness of the model.

Zoom In Zoom Out Reset image size

The trajectories in Figure 3 show that the eccentricity, although initially set to zero, does not remain zero but oscillates about an equilibrium value (Figure 4). We can also adjust the Hamiltonian to obtain greater variation as we move away from the fixed point. Thus, by adjusting the energy of the system through parameters such as semimajor axis, we can control the amplitude of libration. This has implications for tidal heating as discussed in the next section.

Zoom In Zoom Out Reset image size

In addition to amplitude, the period of libration can be chosen based on the relation given in Equation (9). With a set planet mass and set distances such that a 2:1 resonance is achieved, the oscillation period scales with ${m}_{2}^{-2/3}$. We can check the validity of this power law by comparing the libration times we get for various outer moon sizes through simulation. We choose a planet mass of 1.37 M_⊕ and an inner moon eccentricity of 0.094, and plot the time given by Equation (9) alongside the simulation results for five different perturbing masses (Figure 5). The agreement is not exact, but the power law of each differs by less than 10% (0.67 for theory, 0.73 for simulation). The difference can be attributed to the changing eccentricity in the simulation due to being near resonance, whereas the eccentricity was fixed for the theoretical result.

Zoom In Zoom Out Reset image size

The constant phase lag (CPL) model describes the response of a body to a gravitational tug based on the angle, or phase, between the tidal bulge and the line joining the centers of the bodies (e.g., Efroimsky & Makarov 2013). There is a delay between where the planet exerts a gravitational force and where the tidal bulge on the moon appears, and the phase difference is inversely proportional to the quality factor Q that is related to how strongly the body dissipates orbital energy via tides. Although in general the tidal response of a moon can be more complicated than the CPL model, we use this parameterization because it allows us to explore a wide range of system parameters efficiently. For this parameterization, the tidal luminosity is given by

$$\begin{eqnarray}&&{L}_{\mathrm{tide}}=\displaystyle \frac{21{k}_{2}{G}^{3/2}{m}_{p}^{5/2}{R}_{1}^{5}{e}_{1}^{2}}{2{Q}_{1}{a}_{1}^{15/2}},\end{eqnarray} \tag{ 10 }$$

where k₂ is the Love number, R₁ is the radius, and Q₁ is the quality factor of the inner moon, and where we neglect terms of higher order in eccentricity (see Jackson et al. 2008, Henning et al. 2009, and Peters & Turner 2013). Here, we also assume that the moon is tidally locked to its parent planet (Murray & Dermott 1999). Taking the tidal flux to be the luminosity divided by the surface area of the moon and then using the Stefan–Boltzmann law, we obtain for the instantaneous effective temperature of the moon,

$$\begin{eqnarray}&&{T}_{\mathrm{tide}}={\left(\frac{L}{4\pi {\sigma }_{\mathrm{SB}}{R}_{1}^{2}}\right)}^{1/4},\end{eqnarray} \tag{ 11 }$$

where σ_SB is the Stefan–Boltzmann constant.

Returning to the case shown in Figure 4, we simulate the same system in REBOUND and calculate the luminosity and temperature from tides using Equations (10) and (11) and present the results in Figure 6. Because the timescale for the heat transport across the moon is much longer than the period of librations, the surface temperature of the moon will reflect the average rather than the instantaneous heat flux. Taking

$$\begin{eqnarray}&&{L}_{\mathrm{avg}}=\displaystyle \frac{{\int }_{0}^{{t}_{\mathrm{lib}}}{L}_{\mathrm{tide}}\,{dt}}{{t}_{\mathrm{lib}}},\end{eqnarray} \tag{ 12 }$$

we then solve for

$$\begin{eqnarray}&&{T}_{\mathrm{avg}}={\left(\displaystyle \frac{{L}_{\mathrm{avg}}}{4\pi {\sigma }_{\mathrm{SB}}{R}_{1}^{2}}\right)}^{1/4}.\end{eqnarray} \tag{ 13 }$$

L_avg and T_avg are shown by the horizontal red and orange lines, respectively, with T_avg displayed in the top-left corner.

Zoom In Zoom Out Reset image size

The eccentricity at an exact resonance can be analytically estimated for a fixed mass ratio and semimajor axis ratio using Equations (6) and (7) for a constant Hamiltonian K = 0. Solving Equation (6) for x with K and Δ both zero, we get x = 2^1/3. Then, plugging this into Equation (7) with r = 0, a canonical momentum of R = 2^−1/3 is obtained. Inserting this into Equation (4) and solving for the eccentricity, we derive

$$\begin{eqnarray}&&{e}_{1,\mathrm{res}}=0.8{\left[{\left(\displaystyle \frac{{m}_{p}}{{m}_{2}}\right)}^{2/3}-0.76\right]}^{-1/2},\end{eqnarray} \tag{ 14 }$$

and, for m_p ≫ m₂,

$$\begin{eqnarray}&&{e}_{1,\mathrm{res}}\approx 0.1{\left(\displaystyle \frac{{m}_{2}}{{M}_{\mathrm{Io}}}\right)}^{1/3}{\left(\displaystyle \frac{{m}_{p}}{8{M}_{\oplus }}\right)}^{-1/3}.\end{eqnarray} \tag{ 15 }$$

Then, from the eccentricity result, the average tidal heating luminosity temperature is estimated to be

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{res}} & \approx & 2\times {10}^{16}{\left(\displaystyle \frac{Q/{k}_{2}}{{10}^{2}}\right)}^{-1}{\left(\displaystyle \frac{{R}_{1}}{{R}_{\mathrm{Io}}}\right)}^{5}{\left(\displaystyle \frac{{P}_{1}}{2\,\mathrm{day}}\right)}^{-5}\\ & & \times \,{\left(\displaystyle \frac{{m}_{2}}{{M}_{\mathrm{Io}}}\right)}^{2/3}{\left(\displaystyle \frac{{m}_{p}}{8\,{M}_{\oplus }}\right)}^{-2/3}{\rm{W}}.\end{array}\end{eqnarray} \tag{ 16 }$$

Finally, solving for the average temperature due to tides gives

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{res}} & \approx & 303{\left(\displaystyle \frac{Q/{k}_{2}}{{10}^{2}}\right)}^{-1/4}{\left(\displaystyle \frac{{R}_{1}}{{R}_{\mathrm{Io}}}\right)}^{3/4}{\left(\displaystyle \frac{{P}_{1}}{2\,\mathrm{day}}\right)}^{-5/4}\\ & & \times \,{\left(\displaystyle \frac{{m}_{2}}{{M}_{\mathrm{Io}}}\right)}^{1/6}{\left(\displaystyle \frac{{m}_{p}}{8\,{M}_{\oplus }}\right)}^{-1/6}{\rm{K}}.\end{array}\end{eqnarray} \tag{ 17 }$$

Using this method, we estimate the inner moon temperature for the system in Figure 6 to be 369.5 K, within 2% of the simulated result. The difference between the two values can again be attributed to the center of mass location that is offset from the origin. Indeed, this first estimate was for a planet with mass m_p = 4.5 M_⊕ (m_p /m₂ ∼ 300), and increasing the planet-moon mass ratio increases the agreement between the predicted and simulated tidal heating temperatures. Although the constant Hamiltonian configuration described by Equation (2) does not include tidal heating, the formulation described above provides an accurate picture of the short-term heating rates. In Section 3.4, we consider evolution of the system over the longer term through damping timescales and lifetimes of heated moons.

We have also assumed that the inner moon is the test particle in this system as described in Section 2.1 and thus use a massless particle to simulate it in REBOUND. We have analyzed the effect this approximation would have on the dynamics by running similar simulations with an Io-mass moon instead. We find that there is less eccentricity perturbation and thus less tidal heating when using a massive inner moon, an effect that is greater for lower-mass host planets. However, even in the case of a low-mass rocky planet, the difference in tidal heating temperature is less than 20%, and for the Jupiter–Io–Europa system, the difference we find is negligible.

In addition, the CPL tidal heating model is a simplified parameterized method of generalizing tidal interactions based on fixed values Q and k₂. In detail, there is a complicated mechanism of heat transfer within a solid body that depends on parameters such as the viscosity, shear modulus, and density of the material of the mantle, among others. The viscoelastic model of tidal heating aims to incorporate these variables (Dobos et al. 2017). However, this is beyond the scope of this study and will be considered in future work (although, for an example of the impact of using a viscoelastic moon model, see Rovira-Navarro et al. 2021).

Returning to the case of other MMR configurations, we tested other ratios such as 3:1 and 3:2 by running additional simulations and found that a moon in these MMRs will incur a similar amount of heating, depending on the exact orbital periods of the moons. For example, moving the outer moon to a 3 day orbit in the system shown by Figure 6 creates a 3:2 resonance and the inner moon is heated to 338 K. In short, other resonant configurations can be similarly used to trace orbital evolution and tidal heating in exomoons.

Jupiter's moon Io is in a 2:1 resonance with the moon Europa, making it a helpful test case for our heating models. Although the system is in close to perfect resonance, the average eccentricity is significant enough to drive heating within Io. The observed tidal heat flux through Io is 2–4 Wm⁻² (Lainey et al. 2009), corresponding to a surface temperature of 77–92 K (Tyler et al. 2015). To check our heating model against observation, we simulate the three-body system of Jupiter, Io, and Europa and obtain the oscillating tidal luminosity and temperature, which is averaged to obtain a T_avg of 92 K (Figure 7). This agrees with the upper end of observation. Here, we have used k₂ = 0.54 (Lainey et al. 2009) and Q = 36 (Yoder & Peale 1981) for Io. Note the difference in amplitude between the Io case and the zero eccentricity case; the tidal luminosity of Io varies relatively little, and this is because the resonance is close to perfect (the system is near a fixed point).

Zoom In Zoom Out Reset image size

We next explore the range of heating possible in a Jupiter–Io–Europa system to better understand how sensitive the heating is to the Hamiltonian value. We do this by varying initial eccentricity, starting from e = 0.001 and incrementing by 0.001 up to e = 0.008 (Figure 8). The constant Hamiltonian curves encircle the fixed point for all cases. The radius of this circle is largest for e = 0.001 and shrinks as we increase the eccentricity, getting close to the minimum at e = 0.005. As we further increase e, the circle begins to expand again.

Zoom In Zoom Out Reset image size

Figure 8 also includes the heating luminosity and temperature for each value of eccentricity. The trend follows the pattern given by the trajectories of constant Hamiltonian: the variation in heating is a maximum for e = 0.001, decreases as we raise e until e = 0.005, and then increases again. We can tell that the fixed point is between e = 0.001 and e = 0.008 also from the shape of the luminosity curve. At t = 0, the luminosity is at a minimum for e ≤ 0.005 but at a maximum for e ≥ 0.006, signifying that there is minimum variation in heating somewhere between these eccentricities. Although the tidal heating temperature is lowest when we are closest to the fixed point (91 K versus 108 K for e = 0.001), the difference between the maximum and minimum temperatures is less than 20%, displaying that, for a range of Hamiltonian values and even perfect resonance, orbits can give rise to significant tidal heating.

The results above beg the question: What kind of configuration leads to the most or the least heating within a moon? We next explore how central planet mass, and hence composition, and moon orbital period impact tidal heating within an orbiting moon. To study this, we simulate nine different systems with varying m_p and P₁. We choose three different planet masses: 1 M_⊕, 17 M_⊕, and 318 M_⊕, corresponding to Earth, Neptune, and Jupiter. We then choose three pairs of orbital periods for the moons: 2 day/4 day, 4 day/8 day, and 8 day/16 day. The inner moon is massless but has the same radius, Love number, and quality factor as Io, and the perturbing moon is identical except it has the mass of Io. Both moons are in initially circular orbits. Using the same method of calculating tidal luminosity over time, averaging, then finding T_avg, we obtain the results for the nine combinations of these parameters (Figure 9). Here, K and Δ do not take on specific values but are determined by physical parameters as the system evolves. As implied by Equation (16), we get the most heating from low-mass planets and small orbital periods.

Zoom In Zoom Out Reset image size

The maximum amount of heating is for an Earth-like planet with moons in orbits of two and four days. The temperature is 481 K, even greater than the equilibrium temperature of Mercury! This highlights how strongly a moon can be tidally heated when it moves into a resonance—even if it starts in a circular orbit. As the periods are increased across each row, the amount of heating dramatically decreases as expected from the strong scaling with period in Equation (16). Similarly, as planet mass is increased moving down each column, the moon–planet distance again increases to retain the same orbital periods and this drives heating down as well. Another effect evident from Figure 9 is the trend of libration time. As the orbital periods and planet mass increase, the oscillation period for luminosity and temperature increase as predicted by Equation (9).

The large amount of tidal heating found in the previous sections will sap energy from the moon's orbit and eventually move the moon out of resonance. The evolution of the mean motion, n, of the moon depends on the quality factor and Love number of both the planet, Q_p and k_2,p , and the moon. Assuming the planet rotates at a faster rate than the moon orbits, and taking eccentricity to second order, this change in mean motion is given by Boué & Efroimsky (2019):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{n}\displaystyle \frac{{dn}}{{dt}} & = & -\displaystyle \frac{{m}_{1}}{{m}_{p}}{\left(\displaystyle \frac{{R}_{p}}{{a}_{1}}\right)}^{5}\displaystyle \frac{{k}_{2,p}}{{Q}_{p}}\left(\displaystyle \frac{9}{2}+\displaystyle \frac{459}{8}{e}^{2}\right)n\\ & & +\,\displaystyle \frac{171}{2}\displaystyle \frac{{m}_{p}}{{m}_{1}}{\left(\displaystyle \frac{{R}_{1}}{{a}_{1}}\right)}^{5}\displaystyle \frac{{k}_{2}}{Q}{e}^{2}n.\end{array}\end{eqnarray} \tag{ 18 }$$

Similar to Rovira-Navarro et al. (2021), we define a migration timescale of the moon

$$\begin{eqnarray}&&\displaystyle \frac{1}{{t}_{\mathrm{mig}}}=-\displaystyle \frac{9}{2}\displaystyle \frac{{m}_{1}}{{m}_{p}}{\left(\displaystyle \frac{{R}_{p}}{{a}_{1}}\right)}^{5}\displaystyle \frac{{k}_{2,p}}{{Q}_{p}}n,\end{eqnarray} \tag{ 19 }$$

and an eccentricity damping timescale

$$\begin{eqnarray}&&\displaystyle \frac{1}{{t}_{\mathrm{damp}}}=\displaystyle \frac{21}{2}\displaystyle \frac{{m}_{p}}{{m}_{1}}{\left(\displaystyle \frac{{R}_{1}}{{a}_{1}}\right)}^{5}\displaystyle \frac{{k}_{2}}{Q}n,\end{eqnarray} \tag{ 20 }$$

so that Equation (18) becomes

$$\begin{eqnarray}&&\displaystyle \frac{1}{n}\displaystyle \frac{{dn}}{{dt}}=\displaystyle \frac{1}{{t}_{\mathrm{mig}}}\left(1+\displaystyle \frac{51}{4}{e}^{2}\right)+\displaystyle \frac{57}{7}\displaystyle \frac{1}{{t}_{\mathrm{damp}}}{e}^{2}.\end{eqnarray} \tag{ 21 }$$

Thus, when a two-moon system migrates into a 2:1 resonance, we expect the phase of enhanced heating we find to last roughly this damping timescale. Applying Equation (20) to the nine cases presented in Figure 9 provides a sense of the range of damping timescales. At the extremely short-lived end, a moon in a 2 day orbit around a 1 M_⊕ planet peaks at 481 K and will have its eccentricity damped in less than 1 Myr. At the other limit, the eccentricity of a moon in an 8 day orbit around a 1 M_J planet that is heated to only 31 K will last over 100 Myr. For many of our most strongly heated cases, it may be difficult to find the moons in such a state unless a large number of planets are surveyed. That said, finding moons that are strongly heated will put strong constraints on the age and dynamical history of these systems.

Even if the moons move out of resonance, there is still the opportunity for the moons to migrate from tidal interaction with the planet, move back toward resonance, and reinvigorate the heating. To assess the steady-state heating for a moon that experiences this over secular timescales, we follow the arguments of Rovira-Navarro et al. (2021). We summarize the rough order of magnitude conclusions here, but suggest the reader consult this work for a more detailed discussion.

First, for a two-moon system, the inner moon must migrate outward faster than the outer moon to maintain resonance. Comparing these timescales for the two moons and setting 2P₁ = P₂, we find that m₂ ≲ 20m₁ for t_mig,1 < t_mig,2. Second, we estimate the evolution over secular timescales by balancing eccentricity damping (which is trying to move the moons out of resonance) with the outward migration of the inner moon (which is trying to move the moons into resonance). Setting Equation (21) to zero and solving for e gives an equilibrium eccentricity

$$\begin{eqnarray}\begin{array}{rcl}{e}_{\mathrm{eq}} & = & 0.055{\left(\displaystyle \frac{Q}{{Q}_{p}}\right)}^{1/2}{\left(\displaystyle \frac{{k}_{2}}{{k}_{2,p}}\right)}^{-1/2}{\left(\displaystyle \frac{{R}_{1}}{{R}_{\mathrm{Io}}}\right)}^{-5/2}\\ & & \times \,{\left(\displaystyle \frac{{R}_{p}}{2{R}_{\oplus }}\right)}^{5/2}\left(\displaystyle \frac{{m}_{1}}{{M}_{\mathrm{Io}}}\right){\left(\displaystyle \frac{{m}_{p}}{8\,{M}_{\oplus }}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 22 }$$

Substituting this into Equation (10) results in an equilibrium luminosity

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{eq}} & = & {10}^{16}\,\ {\left(\displaystyle \frac{{Q}_{p}/{k}_{2,p}}{100}\right)}^{-1}{\left(\displaystyle \frac{{R}_{p}}{2{R}_{\oplus }}\right)}^{5}\\ & & \times \,{\left(\displaystyle \frac{{P}_{1}}{2\,d}\right)}^{-5}{\left(\displaystyle \frac{{m}_{1}}{{M}_{\mathrm{Io}}}\right)}^{2}{\left(\displaystyle \frac{{m}_{p}}{8\,{M}_{\oplus }}\right)}^{-2}\,{\rm{W}}.\end{array}\end{eqnarray} \tag{ 23 }$$

Such a luminosity could be expected for multimoon systems that are continuously being driven into resonance by migration, representing a longer-term steady-state tidal luminosity.

Low-mass rocky planets generate the most heat within an orbiting moon for a given period ratio (Equations (16) and (17)). However, rocky planets are also the most dissipative, and their moons will undergo the quickest and most significant orbital evolution (Sasaki et al. 2012; Tokadjian & Piro 2020). Furthermore, at least if we use our own solar system as an example, rocky planets may be less likely to harbor compact, multimoon systems. Gas giants such as Jupiter are the least dissipative so that their moons undergo less orbital evolution, but they also cause the least amount of heating. The ideal compromise between moon stability and temperature is thus found in moderate-sized planets such as Neptune that have relatively low dissipation and give rise to appreciable tidal heating. These types of planets can range in size from just larger than super-Earths to ice giants. The former are called sub-Neptunes and have radii of just over 1.7 R_⊕ and masses of typically less than 12 M_⊕ (Lopez & Fortney 2014). These types of planets give the best combination of low mass and slow damping timescale and will be favored moving forward.

Another important consideration is the range of periods where a moon can exist around a given planet. A moon that gets too close to the planet is subject to tidal dissipation and a moon that wanders too far away will be stripped away by the star (e.g., Domingos et al. 2006) and in many cases may eventually collide with planet (Hansen 2022). Starting from the list of all confirmed exoplanets, we categorize those with radius of less than 2.4 R_⊕ or mass of less than 14 M_⊕ as rocky Earth-like planets and those with radius of 1.75–3.5 R_⊕ or mass of 4–14 M_⊕ as Neptune-like, considering those planets that fall under both categories as each case separately. From this list, there are 25 rocky planets with a minimum moon retention time of 300 Myr and 85 Neptune-like planets that can hold a moon for over 5 Gyr. We use the method described in Tokadjian & Piro (2020) to calculate the maximum moon periods for these 110 exoplanets (Figure 10), and set Q_p = 12 for Earth-like planets and Q_p = 10⁴ for Neptune-like planets (Goldreich & Soter 1966; Murray & Dermott 1999). The minimum moon period is set by the radius at which the moon is tidally disrupted by the planet, which depends on the moon's density (Frank et al. 2002). Since we choose an Io-sized moon for all calculations, the value for all cases is 0.24 days. The maximum moon period varies from about 10 days to 90 days, demonstrating that there is a wide range of orbital periods where a moon may exist for all of these planets. Motivated by the Jupiter–Io–Europa system, we choose orbital periods of 2 days and 4 days for moons in resonance for the remainder of this work. This ensures that the inner moon is close enough to the planet to incur significant tidal heating, but far enough to maintain stability in the wake of orbital evolution. Using the scalings we derive analytically and numerically in Section 3, these results can be extrapolated to other moon periods.

Zoom In Zoom Out Reset image size

When trying to identify planets for which exomoon detection is possible, it is also important to consider the contrast between the temperatures of the moon and planet. As we will describe in the next subsection, a larger difference in temperature is helpful in observations. Thus, for each exoplanet candidate, we estimate an irradiation temperature

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}={\left(\displaystyle \frac{{L}_{s}}{16\pi {a}_{p}^{2}{\sigma }_{\mathrm{SB}}}\right)}^{1/4},\end{eqnarray} \tag{ 24 }$$

where L_s is the luminosity of the parent star and a_p is the semimajor axis of the planet. This temperature, along with the planet radius and distance to Earth, are compared for 30 exoplanets in Figure 11. These 30 planets are a sample from the list of planets in Figure 10 made up of the 15 closest Neptune-like planets and 15 closest rocky exoplanets. The vertical lines split the plot into three segments that correspond to planet composition: strictly rocky, rocky or Neptune-like, and strictly Neptune-like based on the radius cut described in Fulton et al. (2017). The horizontal lines represent the habitable zone, while circles and triangles show whether the planet has a determined mass or not, respectively. The color of each marker indicates the distance to Earth.

Zoom In Zoom Out Reset image size

As described above, planets in the sweet spot for moon stability and tidal heating have the radius of a sub-Neptune (toward the right of the plot) and are nearby (blue markers). To maximize planet-moon temperature contrast, planets that have a low mass and radius (toward the left of the plot), are relatively cool in temperature (toward the bottom of the plot), and are nearby (blue markers) would be ideal. For planets with no estimated mass, a mass is assigned based on the radius and density cuts described in Tokadjian & Piro (2020).

Given the relative difficulty of detecting even rocky planets, finding the first rocky exomoon is a challenging task. Thus, novel and indirect methods may be necessary when searching for exomoons. This includes the secondary eclipse method and the search for sodium signatures that would result from volcanic outgassing. In the following, we examine the potential of these methods and follow up with a discussion of promising missions and instruments that may be applicable.

We first study the method of detecting tidally heated exomoons via the technique of secondary transits. This depends on analyzing the relatively small dip in light as the star eclipses the planet. If the planet hosts a moon that is significantly heated, the inferred dayside planet luminosity would be larger than expected.

The amplitude of the secondary eclipse signal is given by the fractional change of light when a planet and moon, with luminosities of L_p and L_m , respectively, go behind the star with luminosity L_s ,

$$\begin{eqnarray}&&{r}_{t}=\displaystyle \frac{{L}_{p}+{L}_{m}}{{L}_{s}}=\displaystyle \frac{{R}_{p}^{2}{T}_{\mathrm{irr}}^{4}+{R}_{1}^{2}({T}_{\mathrm{tide}}^{4}+{T}_{\mathrm{irr}}^{4})}{{R}_{s}^{2}{T}_{s}^{4}},\end{eqnarray} \tag{ 25 }$$

where we take the equilibrium temperature for the planet using Equation (24). Thus, to maximize the signal of the moon in comparison to the planet, it is preferable to consider maximizing the ratios R₁/R_p and T_tide/T_irr. As an example, the blackbody curves in Figure 12 show the specific luminosity of planet and tidally heated inner moon for K2-18b, a 8.88 M_⊕ planet on the cusp between super-Earths and sub-Neptunes. The system has an equilibrium temperature and moon heating temperature of 293 K. It is apparent from this plot that the moon's luminosity is not enough to be detected apart from the planet because here R₁ ≪ R_p and T_tide ≈ T_irr. Most of the sub-Neptunes in our sample face the same issue: because these planets have larger radii than rocky planets and incur relatively little tidal heating in an orbiting moon, the temperature contrast between planet and moon is not significant enough to separate the luminosities of the two bodies.

Zoom In Zoom Out Reset image size

That said, the signal of the moon can be much larger when the parent planet is a smaller sub-Neptune or a rocky Earth-like planet. In particular, we find the light curves of planet and moon begin to separate when the ratio between the moon's peak luminosity and the planet's luminosity at that wavelength is greater than a factor of two. Table 1 summarizes this luminosity ratio along with ${T}_{\mathrm{res}}$ (Equation (17)), distance from Earth, and planet type for 15 exoplanets with the highest luminosity ratio of the 30 shown in Figure 11. Also listed is the moon–planet luminosity ratio using L_eq (Equation (23)), which is relevant for longer-term tidal heating. For the secondary eclipse method, we focus on the short-term heating rates where the moon's luminosity may be much greater than the planet's. The planet with the highest moon-to-planet luminosity ratio is TOI-4353.01, for which the moon outshines the planet by two orders of magnitude at the peak wavelength of the heated moon. Figure 13 shows the blackbody plots for this case: the red curve traces the moon's luminosity, which is clearly distinct from the blue curve of the planet. The tidal heating temperature of the moon is 414 K so that the total temperature of the moon, irradiation plus tidal heating, is 417 K. Using Equation (25) for TOI-4353.01, we get r_t = 0.026 ppm.

Zoom In Zoom Out Reset image size

Table 1. Fifteen Exoplanets with Largest L₁/L_p

Planet	${T}_{\mathrm{res}}$	λ1,peak	${L}_{1,\mathrm{res}}/{l}_{p}$	L1,eq/Lp	Distance	Type
	(K)	(μm)			(pc)
TOI-4353.01	414	6.96	109	0.61	25.1	Rocky
Wolf-1061 d	338	8.54	54	0.63	4.31	Rocky
Wolf-1061 d	338	8.54	24	0.17	4.31	Neptune
TOI-4191.01	320	9.01	18	0.37	83.9	Rocky
TOI-4328.01	354	8.10	7.9	0.18	25.0	Rocky
Kepler-1630 b	326	8.73	5.6	0.04	331	Rocky
GJ 3138 d	320	8.94	4.8	0.16	28.5	Rocky
Kepler-441 b	377	7.56	3.9	0.11	268	Rocky
Kepler-62 f	401	7.07	2.6	0.08	301	Rocky
Kapteyn c	343	8.34	2.3	0.04	3.93	Neptune
GJ 3138 d	320	8.94	2.0	0.16	28.5	Neptune
TOI-4349.01	315	9.00	0.88	0.06	72.5	Rocky
Kepler-1536 b	333	8.53	0.56	0.02	387	Neptune
TOI-4503.01	333	8.26	0.25	0.03	62.5	Rocky
Kepler-452 b	378	7.14	0.25	0.03	552	Rocky

Download table as:  ASCIITypeset image

However, this signal can be made even larger by focusing on a narrower band of wavelengths. By integrating over Planck's Law and normalizing by the disk area, similar to the method described in Charbonneau et al. (2005), the precision for luminosity over the wavelength range λ₁ to λ₂ becomes

$$\begin{eqnarray}&&{r}_{\lambda }=\displaystyle \frac{{\int }_{{\lambda }_{1}}^{{\lambda }_{2}}\left[{R}_{p}^{2}{B}_{\lambda }({T}_{\mathrm{irr}})+{R}_{1}^{2}\left({B}_{\lambda }({T}_{\mathrm{tide}})+{B}_{\lambda }({T}_{\mathrm{irr}})\right)\right]\,d\lambda }{{R}_{s}^{2}{\int }_{{\lambda }_{1}}^{{\lambda }_{2}}{B}_{\lambda }({T}_{s})\,d\lambda }.\end{eqnarray} \tag{ 26 }$$

Here, B_λ (T) is the spectral radiance given by Planck's Law integrated over half the hemisphere to give the specific luminosity. In the case of TOI-4353.01, the moon's radiation peaks at 6.96 μm, so observations in the range of 6–8 μm would be ideal to collect the most light from the moon as possible (Figure 14). Applying Equation (26), we obtain r_λ = 0.419 ppm.

Zoom In Zoom Out Reset image size

Volcanic vents on a tidally heated moon will eject certain kinds of elements such as sodium and potassium, whose clouds may be detectable via transit spectroscopy. In fact, Io is currently losing mass at a rate of about ${\dot{M}}_{\mathrm{Io}}=1000\,\mathrm{kg}$ s⁻¹, a fraction of which is neutral sodium that has been shown to be detectable through sodium line spectroscopy (Oza et al. 2019). Similarly, observations of exoplanets with a tidally heated moon that is venting sodium at a similar or higher rate may provide evidence of the exomoon. Following the methods detailed in Oza et al. (2019), the average column density of sodium viewed against the background star is

$$\begin{eqnarray}&&{N}_{\mathrm{Na}}=\displaystyle \frac{\dot{{M}_{1}}\tau }{{m}_{\mathrm{Na}}\pi {R}_{s}^{2}},\end{eqnarray} \tag{ 27 }$$

where m_Na is the mass of a sodium atom, R_s is the stellar radius, τ is the Na i lifetime, and $\dot{{M}_{1}}$ is the mass-loss rate of sodium from the inner moon due to volcanism. We expect this mass-loss rate to roughly scale proportionally to the energy input from tidal heating. Thus, we scale to Io, giving

$$\begin{eqnarray}&&\dot{{M}_{1}}=\displaystyle \frac{{L}_{1}}{{L}_{\mathrm{Io}}}\dot{{M}_{\mathrm{Io}}}.\end{eqnarray} \tag{ 28 }$$

Here, L₁ and L_Io are the tidal luminosities of the inner moon around the exoplanet and Io, respectively, where the latter is 1.25 × 10¹⁴ W (Veeder et al. 2012). When we evaluate the damping rates for hypothetical moons, we find that they are very short, ∼1 Myr, for some of the most strongly heated cases. This means it may be difficult to build up a sufficiently large cloud of volcanic material and it is unlikely to be present long enough to have a chance of detection. For this reason, we focus on cases where resonance is maintained by migration and use L_eq, given by Equation (23), for estimating L, where Q_p /k_2,p ≈ 300 for Earth-like planets and Q_p /k_2,p ≈ 3000 for Neptune-like planets. Table 2 summarizes the 15 systems with the highest mass-loss rate from the list of nearest and most stable exoplanets. We indicate $\dot{{M}_{1}}$ in both kg s⁻¹ and relative to Io and find that many of these tidally heated moons in resonance are likely to vent sodium at a rate that is tens of times that of Io.

Table 2. Fifteen Exoplanets with Largest Mass-loss Rate

Planet	$\dot{{M}_{1}}$ (kg s−1)	$\dot{{M}_{1}}/\dot{{M}_{\mathrm{Io}}}$	Distance (pc)	Type
TOI-4353.01	1.18 × 105	118	25.1	Rocky
Kepler-62 f	1.11 × 105	111	301	Rocky
Kepler-452 b	0.98 × 105	98.4	552	Rocky
Kepler-441 b	0.98 × 105	98.1	268	Rocky
TOI-4328.01	0.86 × 105	86.1	25	Rocky
TOI-2091.01	0.78 × 105	78.6	70.2	Rocky
Kepler-1040 b	0.77 × 105	76.8	482	Rocky
TOI-4503.01	0.76 × 105	75.9	62.5	Rocky
Kepler-1540 b	0.60 × 105	60.0	245	Neptune
Kepler-1536 b	0.48 × 105	47.4	387	Neptune
HD 216520 c	0.43 × 105	42.6	19.6	Neptune
Kepler-126 d	0.36 × 105	35.7	237	Neptune
Kapteyn c	0.33 × 105	33.6	3.93	Neptune
Wolf-1061 d	0.26 × 105	25.5	4.31	Neptune
HD 136352 d	0.22 × 105	21.9	14.7	Neptune

Download table as:  ASCIITypeset image

Future missions may be able to detect exomoons through the indirect methods we present but must have the necessary sensitivity and a well-matched wavelength range. For example, the planned satellite mission PLAnetary Transits and Oscillations of stars (PLATO; Heng et al. 2014) will have a photometric precision of 50 ppm hr^1/2, which is too large to detect this secondary eclipse of the system TOI-4353.01 (Marchiori et al. 2019). In addition, PLATO will observe at a wavelength range of 0.5–1 μm, and r_λ is much too small for this range to be detectable. However, this wavelength range includes the strong sodium doublet at 589 nm so that sodium line spectroscopy with PLATO would be possible.

The James Webb Space Telescope (JWST) is able to perform spectroscopy using a mid-infrared instrument for the wavelength range 0.6–29 μm, which is just outside the sodium doublet wavelength but fits the ideal case for the secondary eclipse of TOI-4353.01 (Gardner et al. 2006). However, with a median precision of 170 ppm for wavelengths of less than 10 μm, JWST is not as precise as PLATO (Venot et al. 2020) and would not be able to detect the secondary eclipse of TOI-4353.01. Varying potential exomoon size, resonance order, and tidal heating models should be explored in future work to better estimate the prospects for this telescope to find an exomoon.

Similarly, the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) will be sensitive to infrared radiation between 2 and 7.8 μm, the edge of which corresponds to the peak of the test moon around TOI-4353.01 (Pascale et al. 2018). With an estimated photometric precision of 10-100 ppm, ARIEL will have a sensitivity comparable to that of PLATO, and thus might be the best option to date in terms of detecting a tidally heated exomoon around an exoplanet using secondary transits.

Other space telescopes such as the CHaracterising ExOPlanets Satellite (CHEOPS) and the planned Roman Space Telescope (formerly WFIRST) have similar wavelength ranges as PLATO and likely will not find exomoons using the technique of secondary eclipses (Cessa et al. 2017; Johnson et al. 2020). Nevertheless, these missions are set to answer critical questions about extraterrestrial habitability with techniques such as microlensing to discover the smallest planets to date. Smaller planets would incur more heating in their moons for a given resonance mode, but these moons will likely be smaller in size and thus less luminous. Therefore, even with an ideal case of a significantly heated Io-sized moon, it would be difficult to reach the precision needed to detect such a moon with the secondary eclipse method.

Future ground-based observatories such as the Giant Magellan Telescope (GMT), Extremely Large Telescope (ELT), and Thirty Meter Telescope (TMT) will all have instruments that can be used in a wide range of infrared astronomy (DePoy et al. 2012; Thatte et al. 2014; Skidmore 2015). Given that the atmosphere's infrared window is at slightly longer wavelengths (∼8–14 μm versus a peak of 7.35 μm), the exomoon in question would need to have a lower overall temperature. This is more likely for larger exoplanets, such as the sub-Neptune class mentioned earlier, so long as the equilibrium temperature of the bodies is small enough for the moon to still outshine its parent planet. Because these observatories will also observe in the optical, sodium line spectroscopy of potential tidally heated moons is promising.